// Gauss-Kronrod-Patterson quadrature coefficients for use in
// quadpack routine qng. These coefficients were calculated with
// 101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov
// 1981. 

namespace integ {

    static const int NGKPLEVELS = 7;

    inline int gkp_n(int level) 
    { 
        assert(level >= 0 && level < NGKPLEVELS);
        static const int ngkp[NGKPLEVELS] = {0,1,2,4,8,16,32};
        return ngkp[level]; 
    }

    template <class T> 
    inline const std::vector<T>& gkp_x(int level)
    {
        static const std::vector<T> vx1(0);

        static const T ax3[1] = {
            0.7745966692414833770358531
        };
        static const std::vector<T> vx3(ax3,ax3+1);

        static const T ax7[2] = {
            0.9604912687080202834235071,
            0.4342437493468025580020715
        };
        static const std::vector<T> vx7(ax7,ax7+2);

        static const T ax15[4] = {
            0.9938319632127550222085128,
            0.8884592328722569988904202,
            0.6211029467372264029406874,
            0.2233866864289668816282040
        };
        static const std::vector<T> vx15(ax15,ax15+4);

        static const T ax31[8] = {
            0.9990981249676675976622261,
            0.9815311495537401068673619,
            0.9296548574297400566701257,
            0.8367259381688687355027538,
            0.7024962064915270786098002,
            0.5313197436443756239721034,
            0.3311353932579768330926408,
            0.1124889431331866257458433
        };
        static const std::vector<T> vx31(ax31,ax31+8);

        static const T ax63[16] = {
            0.9998728881203576119379568,
            0.9972062593722219590764525,
            0.9886847575474294799385289,
            0.9721828747485817965780588,
            0.9463428583734029051484962,
            0.9103711569570042924977907,
            0.8639079381936904771464159,
            0.8069405319502176118563080,
            0.7397560443526947586772178,
            0.6629096600247805954610153,
            0.5771957100520458148436910,
            0.4836180269458410275621533,
            0.3833593241987303469164852,
            0.2777498220218243150653564,
            0.1682352515522074649823133,
            0.05634431304659278997196786
        };
        static const std::vector<T> vx63(ax63,ax63+16);

        static const T ax127[32] = {
            0.9999824303548916002883871,
            0.9995987996719106826208247,
            0.9983166353184073925744239,
            0.9957241046984071885071788,
            0.9914957211781061323986149,
            0.9853714995985203711137519,
            0.9771415146397057141563962,
            0.9666378515584165670922798,
            0.9537300064257611364147486,
            0.9383203977795928836548223,
            0.9203400254700124207298214,
            0.8997448997769400366386332,
            0.8765134144847052697416266,
            0.8506444947683502797578274,
            0.8221562543649804073725271,
            0.7910849337998483614346381,
            0.7574839663805136379262696,
            0.7214230853700989154849762,
            0.6829874310910792280870776,
            0.6422766425097595137741136,
            0.5994039302422428929742510,
            0.5544951326319325488663814,
            0.5076877575337166021547831,
            0.4591300119898323328735020,
            0.4089798212298886724090317,
            0.3574038378315321523762149,
            0.3045764415567140433353240,
            0.2506787303034831766129571,
            0.1958975027111001539154602,
            0.1404242331525601745938196,
            0.08445404008371088371018217,
            0.02818464894974569433939733
        };
        static const std::vector<T> vx127(ax127,ax127+32);

        static const std::vector<T>* x[NGKPLEVELS] = 
        { &vx1,&vx3,&vx7,&vx15,&vx31,&vx63,&vx127 };

        assert(level >= 0 && level < NGKPLEVELS);
        return *x[level];
    }

    inline const std::vector<T>& gkp_wa(int level)
    {
        static const std::vector<T> vw3a(0);

        static const T aw7a[1] = {
            0.2684880898683334407285693
        };
        static const std::vector<T> vw7a(aw7a,aw7a+1);

        static const T aw15a[3] = {
            0.1344152552437842203599688,
            0.05160328299707973969692012,
            0.2006285293769890210339319
        };
        static const std::vector<T> vw15a(aw15a,aw15a+3);

        static const T aw31a[7] = {
            0.06720775429599070354040106,
            0.02580759809617665356464612,
            0.1003142786117955787712936,
            0.008434565739321106246314930,
            0.04646289326175798654140464,
            0.08575592004999035115418652,
            0.1095784210559246382366884
        };
        static const std::vector<T> vw31a(aw31a,aw31a+7);

        static const T aw63a[15] = {
            0.03360387714820773054173399,
            0.01290380010035126562597665,
            0.05015713930589953741367955,
            0.004217630441558854839084227,
            0.02323144663991026944325649,
            0.04287796002500773449291230,
            0.05478921052796286503221753,
            0.001265156556230068011372609,
            0.008223007957235929669257784,
            0.01797855156812827033289605,
            0.02848975474583354861250609,
            0.03843981024945553203864035,
            0.04681355499062801240264808,
            0.05283494679011651986207666,
            0.05597843651047631940755338
        };
        static const std::vector<T> vw63a(aw63a,aw63a+15);

        static const T aw127a[31] = {
            0.01680193857410386527086942,
            0.006451900050175736922805087,
            0.02507856965294976870683977,
            0.002108815245726632878535242,
            0.01161572331995513472698495,
            0.02143898001250386724645616,
            0.02739460526398143251610877,
            0.0006326073193626332640525769,
            0.004111503978654693047167841,
            0.008989275784064135723280604,
            0.01424487737291677430634157,
            0.01921990512472776601932028,
            0.02340677749531400620132404,
            0.02641747339505825993103833,
            0.02798921825523815970377669,
            0.0001807395644453909103150119,
            0.001289524082610417407932182,
            0.003057753410175531136173120,
            0.005249123454808859125133999,
            0.007703375233279741848165979,
            0.01029711695795635552368646,
            0.01293483966360737345473396,
            0.01553677555584398243992842,
            0.01803221639039128632005310,
            0.02035775505847215946694702,
            0.02245726582681609870712712,
            0.02428216520333659935797356,
            0.02579162697602422938840455,
            0.02695274966763303196343848,
            0.02774070217827968199391920,
            0.02813884991562715063629767
        };
        static const std::vector<T> vw127a(aw127a,aw127a+31);

        static const std::vector<T>* wa[NGKPLEVELS] = 
        {0,&vw3a,&vw7a,&vw15a,&vw31a,&vw63a,&vw127a};

        assert(level >= 1 && level < NGKPLEVELS);
        return *wa[level];
    }

    inline const std::vector<T>& gkp_wb(int level)
    {
        static const T aw1b[1] = {
            2.000000000000000000000000
        };
        static const std::vector<T> vw1b(aw1b,aw1b+1);

        static const T aw3b[2] = {
            0.5555555555555555555555556,
            0.8888888888888888888888889
        };
        static const std::vector<T> vw3b(aw3b,aw3b+2);

        static const T aw7b[3] = {
            0.1046562260264672651938239,
            0.4013974147759622229050518,
            0.4509165386584741423451101
        };
        static const std::vector<T> vw7b(aw7b,aw7b+3);

        static const T aw15b[5] = {
            0.01700171962994026033902742,
            0.09292719531512453768589422,
            0.1715119091363913807873532,
            0.2191568584015874964036932,
            0.2255104997982066873864225
        };
        static const std::vector<T> vw15b(aw15b,aw15b+5);

        static const T aw31b[9] = {
            0.002544780791561874415402782,
            0.01644604985438781093378839,
            0.03595710330712932209677783,
            0.05697950949412335741219737,
            0.07687962049900353104270519,
            0.09362710998126447361665878,
            0.1056698935802348097438159,
            0.1119568730209534568801436,
            0.1127552567207686916071499
        };
        static const std::vector<T> vw31b(aw31b,aw31b+9);

        static const T aw63b[17] = {
            0.0003632214818455306596935806,
            0.002579049794685688272427796,
            0.006115506822117246339678284,
            0.01049824690962132189827284,
            0.01540675046655949780213083,
            0.02059423391591271114918856,
            0.02586967932721474691075827,
            0.03107355111168796487988439,
            0.03606443278078257264010716,
            0.04071551011694431893389410,
            0.04491453165363219741425425,
            0.04856433040667319871594712,
            0.05158325395204845877680910,
            0.05390549933526606392687695,
            0.05548140435655936398783841,
            0.05627769983125430127259535,
            0.05637762836038471738766256
        };
        static const std::vector<T> vw63b(aw63b,aw63b+17);

        static const T aw127b[33] = {
            0.00005053609520786089516638322,
            0.0003777466463269839352971478,
            0.0009383698485423815640378046,
            0.001681142865421469902932108,
            0.002568764943794020373297319,
            0.003572892783517299649375536,
            0.004671050372114321747405966,
            0.005843449875835639507559476,
            0.007072489995433555468046319,
            0.008342838753968157705584124,
            0.009641177729702536695298303,
            0.01095573338783790164803273,
            0.01227583056008277008696633,
            0.01359157100976554678957292,
            0.01489364166481518203481040,
            0.01617321872957771994194796,
            0.01742193015946417374715226,
            0.01863184825613879018631404,
            0.01979549504809749948802772,
            0.02090585144581202385222185,
            0.02195636630531782493926050,
            0.02294096422938774876080053,
            0.02385405210603854008044603,
            0.02469052474448767690906084,
            0.02544576996546476581257440,
            0.02611567337670609768049881,
            0.02669662292745035990615470,
            0.02718551322962479181920860,
            0.02757974956648187303486871,
            0.02787725147661370160852380,
            0.02807645579381724660684785,
            0.02817631903301660213065358,
            0.02818881418019235869383128
        };
        static const std::vector<T> vw127b(aw127b,aw127b+33);

        static const std::vector<T>* wb[NGKPLEVELS] = 
        { &vw1b,&vw3b,&vw7b,&vw15b,&vw31b,&vw63b,&vw127b };

        assert(level >= 0 && level < NGKPLEVELS);
        return *wb[level];
    }

#ifdef EXTRA_PREC_H
    template <> 
    inline const std::vector<Quad>& gkp_x<Quad>(int level)
    {
        static const std::vector<Quad> vx1(0);

        static const Quad ax3[1] = {
            "0.7745966692414833770358531"
        };
        static const std::vector<Quad> vx3(ax3,ax3+1);

        static const Quad ax7[2] = {
            "0.9604912687080202834235071",
            "0.4342437493468025580020715"
        };
        static const std::vector<Quad> vx7(ax7,ax7+2);

        static const Quad ax15[4] = {
            "0.9938319632127550222085128",
            "0.8884592328722569988904202",
            "0.6211029467372264029406874",
            "0.2233866864289668816282040"
        };
        static const std::vector<Quad> vx15(ax15,ax15+4);

        static const Quad ax31[8] = {
            "0.9990981249676675976622261",
            "0.9815311495537401068673619",
            "0.9296548574297400566701257",
            "0.8367259381688687355027538",
            "0.7024962064915270786098002",
            "0.5313197436443756239721034",
            "0.3311353932579768330926408",
            "0.1124889431331866257458433"
        };
        static const std::vector<Quad> vx31(ax31,ax31+8);

        static const Quad ax63[16] = {
            "0.9998728881203576119379568",
            "0.9972062593722219590764525",
            "0.9886847575474294799385289",
            "0.9721828747485817965780588",
            "0.9463428583734029051484962",
            "0.9103711569570042924977907",
            "0.8639079381936904771464159",
            "0.8069405319502176118563080",
            "0.7397560443526947586772178",
            "0.6629096600247805954610153",
            "0.5771957100520458148436910",
            "0.4836180269458410275621533",
            "0.3833593241987303469164852",
            "0.2777498220218243150653564",
            "0.1682352515522074649823133",
            "0.05634431304659278997196786"
        };
        static const std::vector<Quad> vx63(ax63,ax63+16);

        static const Quad ax127[32] = {
            "0.9999824303548916002883871",
            "0.9995987996719106826208247",
            "0.9983166353184073925744239",
            "0.9957241046984071885071788",
            "0.9914957211781061323986149",
            "0.9853714995985203711137519",
            "0.9771415146397057141563962",
            "0.9666378515584165670922798",
            "0.9537300064257611364147486",
            "0.9383203977795928836548223",
            "0.9203400254700124207298214",
            "0.8997448997769400366386332",
            "0.8765134144847052697416266",
            "0.8506444947683502797578274",
            "0.8221562543649804073725271",
            "0.7910849337998483614346381",
            "0.7574839663805136379262696",
            "0.7214230853700989154849762",
            "0.6829874310910792280870776",
            "0.6422766425097595137741136",
            "0.5994039302422428929742510",
            "0.5544951326319325488663814",
            "0.5076877575337166021547831",
            "0.4591300119898323328735020",
            "0.4089798212298886724090317",
            "0.3574038378315321523762149",
            "0.3045764415567140433353240",
            "0.2506787303034831766129571",
            "0.1958975027111001539154602",
            "0.1404242331525601745938196",
            "0.08445404008371088371018217",
            "0.02818464894974569433939733"
        };
        static const std::vector<Quad> vx127(ax127,ax127+32);

        static const std::vector<Quad>* x[NGKPLEVELS] = 
        { &vx1,&vx3,&vx7,&vx15,&vx31,&vx63,&vx127 };

        assert(level >= 0 && level < NGKPLEVELS);
        return *x[level];
    }

    template <> 
    inline const std::vector<Quad>& gkp_wa<Quad>(int level)
    {
        static const std::vector<Quad> vw3a(0);

        static const Quad aw7a[1] = {
            "0.2684880898683334407285693"
        };
        static const std::vector<Quad> vw7a(aw7a,aw7a+1);

        static const Quad aw15a[3] = {
            "0.1344152552437842203599688",
            "0.05160328299707973969692012",
            "0.2006285293769890210339319"
        };
        static const std::vector<Quad> vw15a(aw15a,aw15a+3);

        static const Quad aw31a[7] = {
            "0.06720775429599070354040106",
            "0.02580759809617665356464612",
            "0.1003142786117955787712936",
            "0.008434565739321106246314930",
            "0.04646289326175798654140464",
            "0.08575592004999035115418652",
            "0.1095784210559246382366884"
        };
        static const std::vector<Quad> vw31a(aw31a,aw31a+7);

        static const Quad aw63a[15] = {
            "0.03360387714820773054173399",
            "0.01290380010035126562597665",
            "0.05015713930589953741367955",
            "0.004217630441558854839084227",
            "0.02323144663991026944325649",
            "0.04287796002500773449291230",
            "0.05478921052796286503221753",
            "0.001265156556230068011372609",
            "0.008223007957235929669257784",
            "0.01797855156812827033289605",
            "0.02848975474583354861250609",
            "0.03843981024945553203864035",
            "0.04681355499062801240264808",
            "0.05283494679011651986207666",
            "0.05597843651047631940755338"
        };
        static const std::vector<Quad> vw63a(aw63a,aw63a+15);

        static const Quad aw127a[31] = {
            "0.01680193857410386527086942",
            "0.006451900050175736922805087",
            "0.02507856965294976870683977",
            "0.002108815245726632878535242",
            "0.01161572331995513472698495",
            "0.02143898001250386724645616",
            "0.02739460526398143251610877",
            "0.0006326073193626332640525769",
            "0.004111503978654693047167841",
            "0.008989275784064135723280604",
            "0.01424487737291677430634157",
            "0.01921990512472776601932028",
            "0.02340677749531400620132404",
            "0.02641747339505825993103833",
            "0.02798921825523815970377669",
            "0.0001807395644453909103150119",
            "0.001289524082610417407932182",
            "0.003057753410175531136173120",
            "0.005249123454808859125133999",
            "0.007703375233279741848165979",
            "0.01029711695795635552368646",
            "0.01293483966360737345473396",
            "0.01553677555584398243992842",
            "0.01803221639039128632005310",
            "0.02035775505847215946694702",
            "0.02245726582681609870712712",
            "0.02428216520333659935797356",
            "0.02579162697602422938840455",
            "0.02695274966763303196343848",
            "0.02774070217827968199391920",
            "0.02813884991562715063629767"
        };
        static const std::vector<Quad> vw127a(aw127a,aw127a+31);

        static const std::vector<Quad>* wa[NGKPLEVELS] = 
        {0,&vw3a,&vw7a,&vw15a,&vw31a,&vw63a,&vw127a};

        assert(level >= 1 && level < NGKPLEVELS);
        return *wa[level];
    }

    template <> 
    inline const std::vector<Quad>& gkp_wb<Quad>(int level)
    {
        static const Quad aw1b[1] = {
            "2.000000000000000000000000"
        };
        static const std::vector<Quad> vw1b(aw1b,aw1b+1);

        static const Quad aw3b[2] = {
            "0.5555555555555555555555556",
            "0.8888888888888888888888889"
        };
        static const std::vector<Quad> vw3b(aw3b,aw3b+2);

        static const Quad aw7b[3] = {
            "0.1046562260264672651938239",
            "0.4013974147759622229050518",
            "0.4509165386584741423451101"
        };
        static const std::vector<Quad> vw7b(aw7b,aw7b+3);

        static const Quad aw15b[5] = {
            "0.01700171962994026033902742",
            "0.09292719531512453768589422",
            "0.1715119091363913807873532",
            "0.2191568584015874964036932",
            "0.2255104997982066873864225"
        };
        static const std::vector<Quad> vw15b(aw15b,aw15b+5);

        static const Quad aw31b[9] = {
            "0.002544780791561874415402782",
            "0.01644604985438781093378839",
            "0.03595710330712932209677783",
            "0.05697950949412335741219737",
            "0.07687962049900353104270519",
            "0.09362710998126447361665878",
            "0.1056698935802348097438159",
            "0.1119568730209534568801436",
            "0.1127552567207686916071499"
        };
        static const std::vector<Quad> vw31b(aw31b,aw31b+9);

        static const Quad aw63b[17] = {
            "0.0003632214818455306596935806",
            "0.002579049794685688272427796",
            "0.006115506822117246339678284",
            "0.01049824690962132189827284",
            "0.01540675046655949780213083",
            "0.02059423391591271114918856",
            "0.02586967932721474691075827",
            "0.03107355111168796487988439",
            "0.03606443278078257264010716",
            "0.04071551011694431893389410",
            "0.04491453165363219741425425",
            "0.04856433040667319871594712",
            "0.05158325395204845877680910",
            "0.05390549933526606392687695",
            "0.05548140435655936398783841",
            "0.05627769983125430127259535",
            "0.05637762836038471738766256"
        };
        static const std::vector<Quad> vw63b(aw63b,aw63b+17);

        static const Quad aw127b[33] = {
            "0.00005053609520786089516638322",
            "0.0003777466463269839352971478",
            "0.0009383698485423815640378046",
            "0.001681142865421469902932108",
            "0.002568764943794020373297319",
            "0.003572892783517299649375536",
            "0.004671050372114321747405966",
            "0.005843449875835639507559476",
            "0.007072489995433555468046319",
            "0.008342838753968157705584124",
            "0.009641177729702536695298303",
            "0.01095573338783790164803273",
            "0.01227583056008277008696633",
            "0.01359157100976554678957292",
            "0.01489364166481518203481040",
            "0.01617321872957771994194796",
            "0.01742193015946417374715226",
            "0.01863184825613879018631404",
            "0.01979549504809749948802772",
            "0.02090585144581202385222185",
            "0.02195636630531782493926050",
            "0.02294096422938774876080053",
            "0.02385405210603854008044603",
            "0.02469052474448767690906084",
            "0.02544576996546476581257440",
            "0.02611567337670609768049881",
            "0.02669662292745035990615470",
            "0.02718551322962479181920860",
            "0.02757974956648187303486871",
            "0.02787725147661370160852380",
            "0.02807645579381724660684785",
            "0.02817631903301660213065358",
            "0.02818881418019235869383128"
        };
        static const std::vector<Quad> vw127b(aw127b,aw127b+33);

        static const std::vector<Quad>* wb[NGKPLEVELS] = 
        { &vw1b,&vw3b,&vw7b,&vw15b,&vw31b,&vw63b,&vw127b };

        assert(level >= 0 && level < NGKPLEVELS);
        return *wb[level];
    }
#endif

}
